#include "gmshImport.hpp"
#include "meshFunctions.hpp"
#include "topology.hpp"
#include "geometry.hpp"
#include "util.hpp"
#include "streamAccess.hpp"
#include "exceptions.hpp"
#include "timer.hpp"
#include "macros.hpp"
#include "customReductions.hpp"

#include <iomanip>
#include <iterator>
#include <iostream>
#include <unordered_map>
#include <tuple>

namespace hms::gmsh
{

Mesh gmshToMesh( const fs::path& file ){
	std::ifstream is { getReadFile(file, false) };
	return gmshToHms(is);
}

/**
 * @brief Imports a MSH file generated by gmsh into the hms mesh format.
 * Boundary names must be given in the .geo file as:
 * Physical Curve("yourBoundaryName") = { <curve tags (int)> };
 * ONLY boundary patches can have physical names, otherwise they are mistakenly
 * identified as boundary patches!
 * When physical names are stated, gmsh defaults to not saving anything that 
 * is not grouped into a physical entity. To avoid this, please put the 
 * following line into your .geo file:
 * Mesh.SaveAll=1;
 * Alternatively, you may also use the command line option -save_all for 
 * generating the MSH file. For further information, please consult the gmsh
 * documentation, specifically 
 * http://gmsh.info/doc/texinfo/gmsh.html#Elementary-entities-vs-physical-groups
 * @param file the MSH file as ifstream
 * @return a fully initialised mesh
 */
/* notes:
 * No duplicate nodes are stored. They are stored under the lowest EntityDim
 * which uses them. To use them for the mesh, these specifics can be ignored.
 * Instead, only their coordinates and tags are aggregated and stored.
 * Elements: triangles and quadrilaterals form cells. No lines for non-boundary
 * edges are stored in the MSH file, so they need to be generated. Boundary
 * edges are stored as lines, but other lines may be stored as well, so the
 * physical names section should be consulted to see which entities belong to
 * boundary patches and in turn, which elements belong to those entities.
 * These entities themselves are unimportant, because point entities do not need
 * representation outside of being stored as a node already, curves only matter
 * when they belong to a boundary patch, and as there may be multiple curves in
 * one boundary patch, the single curve is eventually unimportant. And there's 
 * only one surface, so no distinction is necessary or would even be useful if
 * there were more.
 * */
Mesh gmshToHms( std::istream& file ){
	const MSH msh { readFile( file ) };
	if ( !msh.isValid )
		throw MeshError( "(gmsh import) Invalid or empty MSH file!");
	
	/* the file is no longer needed */
	auto f { dynamic_cast<std::ifstream*>(&file) };
	/* make sure the passed in istream was actually a file stream */
	if (f) /* and if so, close the file */
		f->close();
	
	/* everything between START and STOP can be executed out of order */
	/* START */
	/**************************************************************************/
	/* collect all node tags associated with cells */
	/**************************************************************************/
	Index nDomainCells { gmsh::nDomainCells(msh) };
	
	Array4Xi      domainCellVertices;
	ArrayX<char> nDomainCellVertices;
	
	std::tie(domainCellVertices, nDomainCellVertices) =
		collectDomainCellNodeTags( msh, nDomainCells );
	
	/**************************************************************************/
	/* collect edges associated with each boundary patch */
	/**************************************************************************/
	/* the boundary edges are created by gmsh, and if physical names are given
	 * to them, the edges can be attributed to their boundary patch. */
	IdxVector boundaryPhysicalTags { boundaryPatchTags(msh) };
	
	Array2Xi boundaryEdges;
	std::vector<ArrayXi> edgeIndicesPerPatch;
	
	std::tie ( boundaryEdges, edgeIndicesPerPatch ) =
		boundaryEdgesPerPatch(msh, boundaryPhysicalTags);
	
	/**************************************************************************/
	/* collect all node coordinates and their tags */
	/**************************************************************************/
	IdxVector nodeTags;
	Matrix2Xs nodes;
	
	std::tie (nodeTags, nodes) = collectNodes(msh);
	
	/* STOP */
	/* everything between START and STOP can be executed out of order */
	
	/**************************************************************************/
	/* translate node tags to indices */
	/**************************************************************************/
	tagToIndex( domainCellVertices, nDomainCellVertices, nodeTags );
	tagToIndex<2>(boundaryEdges, nodeTags);
	
	/**************************************************************************/
	/* create and store domain cell edges */
	/**************************************************************************/
	// hms::Mesh mesh { hms::cellVerticesToMesh(
	// 	std::move(nodes), std::move(nDomainCellVertices), domainCellVertices
	// ) };
	
	// /* now the edge indices for boundary patches must be offset, because they
	//  * occupy the rightmost columns of edgeVertices */
	// for ( Ref<ArrayXi> indices : edgeIndicesPerPatch ){
	// 	indices.array() += mesh.nInnerEdges();
	// }

	Array2Xi edgeVertices;
	Array4Xi domainCellEdges;
	std::tie( edgeVertices, domainCellEdges ) = createDomainCellEdges(
		msh, domainCellVertices, nDomainCellVertices, boundaryEdges
	);
	
	/* convenience */
	Index nEdges         { edgeVertices.cols() };
	Index nBoundaryEdges { boundaryEdges.cols() };
	Index nInnerEdges    { nEdges - nBoundaryEdges };
	
	/* now the edge indices for boundary patches must be offset, because they
	 * occupy the rightmost columns of edgeVertices */
	for ( Ref<ArrayXi> indices : edgeIndicesPerPatch ){
		indices.array() += nInnerEdges;
	}
	
	/**************************************************************************/
	/* store parent cell indices for edges */
	/**************************************************************************/
	Array2Xi adjacentCells { hms::findAdjacentCells(
		nDomainCellVertices, domainCellEdges, edgeVertices.cols()
	) };
	
	/* ghost cells as parents for boundary edges */
	for ( Index i{0}; i<nBoundaryEdges; ++i )
		adjacentCells(1, nInnerEdges + i) = nDomainCells + i;
	
	/**************************************************************************/
	/* create mesh */
	/**************************************************************************/
	hms::Mesh mesh {
		std::move(nodes),
		nInnerEdges,
		std::move(edgeVertices),
		std::move(adjacentCells),
		std::move(nDomainCellVertices),
		std::move(domainCellEdges),
		domainCellVertices
	};
	
	/**************************************************************************/
	/* add boundary patches */
	/**************************************************************************/
	for( uIndex i{0}; i<boundaryPhysicalTags.size(); ++i ){
		for ( const PhysicalName& p : msh.physicalNames ){
			if ( p.dim == Dim::curve && p.tag == boundaryPhysicalTags[i] )
				mesh.addBoundaryPatch(
					p.name, std::move(edgeIndicesPerPatch[i]) );
		}
	}
	
	return mesh;
}



auto nDomainCells( const MSH& msh ) -> Index {
	Index nCells {0};
	OMP_PRAGMA(parallel for reduction (+:nCells) schedule(static))
	for ( const ElementBlock& eb : msh.elementBlocks ){
		if (
			eb.elementType == ElementType::triangle ||
			eb.elementType == ElementType::quadrilateral
		)
			nCells += eb.elementTags.size();
	}
	return nCells;
}



/**
 * @brief The formula for the number of inner edges n_i is
 * n_i = 1/2 * ( 3 * n_t + 4 * n_q - n_b )
 * with n_t and n_q being the number of triangular and quadrilateral cells,
 * respectively, and n_b being the number of boundary edges
 * @param msh
 * @return the number of inner edges
 */
Index nInnerEdges( const MSH& msh, Index nBoundaryEdges ){
	Index nTriangles {0}, nQuadrilaterals {0};
	
	OMP_PRAGMA(
		parallel for reduction (+:nTriangles,nQuadrilaterals) schedule(static)
	)
	for ( const ElementBlock& eb : msh.elementBlocks ){
		if ( eb.elementType == ElementType::triangle )
			nTriangles += eb.elementTags.size();
		if ( eb.elementType == ElementType::quadrilateral )
			nQuadrilaterals += eb.elementTags.size();
	}
	
	return ( 3 * nTriangles + 4 * nQuadrilaterals - nBoundaryEdges ) / 2;
}



auto collectNodes( const MSH& msh )
-> std::tuple<IdxVector, Matrix2Xs>
{
	/* node (i.e. vertex) coordinates and tags */
	Matrix2Xs nodes (2, msh.nNodes);
	IdxVector nodeTags;
	nodeTags.reserve( msh.nNodes );
	
	/* consolidate nodes */
	Index start {0};
	for ( const NodeBlock& nb : msh.nodeBlocks ){
		std::copy(
			nb.nodeTags.begin(), nb.nodeTags.end(), std::back_inserter(nodeTags)
		);
		/* z-coordinate is being ignored */
		nodes.block( 0, start, 2, nb.coords.cols() ) = nb.coords.topRows(2);
		start += nb.coords.cols();
	}
	
	return { nodeTags, nodes };
}



auto collectDomainCellNodeTags( const MSH& msh, Index nDomainCells )
-> std::pair<
	Array4Xi,    // vertex indices
	ArrayX<char> // number of vertices for each cell
>{
	Array4Xi domainCellNodes (4, nDomainCells);
	ArrayX<char>  nDomainCellNodes (   nDomainCells);
	
	Index cell {0};
	for ( const ElementBlock& eb : msh.elementBlocks ){
		if (
			eb.elementType == ElementType::triangle ||
			eb.elementType == ElementType::quadrilateral
		){
			Index nVertices { nNodeTags( eb.elementType ) };
			nDomainCellNodes.segment( cell, eb.nodeTags.cols() )
				= ArrayX<char>::Constant( eb.nodeTags.cols(), nVertices );
			
			/* each column in eb.nodeTags holds the vertices of a cell */
			domainCellNodes.block( 0, cell, nVertices, eb.nodeTags.cols() )
				= eb.nodeTags;
			
			cell += eb.nodeTags.cols();
		}
	}
	return { domainCellNodes, nDomainCellNodes };
}



auto boundaryPatchTags( const MSH& msh ) -> IdxVector {
/* approach for boundary edges:
 * physical names are looked through to find ones that have the dimension of a
 * curve. It is assumed that named curves or curve groups are boundaries.
 * They are stored in a vector, with each index belonging to one boundary patch.
 * Then entity tags of curves belonging to each boundary patch are stored, so
 * that later, elements (i.e. lines) which reference these entities (i.e. 
 * curves) can be connected to that boundary patch.
 * */
	/* look through physical names and store the tags of boundary patch curves
	 * */
	IdxVector boundaryPatchTags;
	for ( const PhysicalName& pn : msh.physicalNames ){
		if ( pn.dim == Dim::curve )
			boundaryPatchTags.push_back( pn.tag );
	}
	return boundaryPatchTags;
}
		


/**
 * @brief Returns node and edge tags for all boundary patches
 * @param msh The representation of a MSH file
 * @param boundaryPatchTags just the tags of those physical names that are
 * boundary patches
 * @return a std::tuple containing three objects: The first is a 2*n matrix
 * holding the node tags for each edge (i.e. tags of their two endpoints). Each
 * column of that matrix represents one each. The other two 
 * 
 * The association between boundary edges and boundary patches/patch
 * names is through 1D-entities, i.e. curves, with a presumed 1-n association.
 * The MSH format allows n-n associations, but those do not make sense for 
 * boundary patches.
 * This routine looks at the association of each curve entity with a physical
 * name, then at all boundary edges and their association with curves, to
 * directly associate edges with boundary patches.
 * Boundary edges cannot be retrieved independently of physical names, because
 * it is possible to create inner edges in gmsh, e.g. to control the resolution
 * along an interface. These would be found as well and mistakenly identified as
 * boundary edges.
 */
auto boundaryEdgesPerPatch( const MSH& msh, const IdxVector& boundaryPatchTags )
-> std::pair< Array2Xi, std::vector<ArrayXi> >
{
	Index nBoundaryPatches { static_cast<Index>( boundaryPatchTags.size() ) };
	/* look through curves and store the entity tags of the curves that belong
	 * to a boundary patch curve at that index */
	std::vector<IdxVector> boundaryCurveTags ( nBoundaryPatches );
	for ( Index i{0}; i<nBoundaryPatches ; ++i ){
		for ( const Curve& c : msh.curves ){
			if ( hms::hasEntry(
					c.physicalTags,
					static_cast<Index>(boundaryPatchTags[i])
				) )
				boundaryCurveTags[i].push_back( c.tag );
		}
	}
	
	/* To associate edges with boundary patches, the number of edges per patch
	 * is calculated as the sum of edges of the curves associated with each
	 * patch */
	ArrayXi nBoundaryEdgesPerPatch { /*sums must be initialised to zero*/
		ArrayXi::Constant( nBoundaryPatches, 0 )
	};
	OMP_PRAGMA(parallel for collapse(2) reduction(+:nBoundaryEdgesPerPatch))
	for ( Index i=0; i<nBoundaryPatches; ++i ){
		for ( const ElementBlock& eb : msh.elementBlocks ){
			if (
				eb.entityDim == Dim::curve &&
				eb.elementType == ElementType::line &&
				hms::hasEntry( boundaryCurveTags[i], eb.entityTag )
			){
				nBoundaryEdgesPerPatch[i] += eb.elementTags.size();
			}
		}
	}
	
	/* The total number of boundaryEdges is now simply the sum of boundary edges
	 * per patch */
	Index nBoundaryEdges { nBoundaryEdgesPerPatch.sum() };
	
	/* look through edges and store the endpoint node tags of those that are
	 * part of boundary patches. By having the outer loop iterating over 
	 * boundary patches, the columns (holding node tags for edge endpoints) are
	 * necessarily ordered by boundary patches, which is nice :)
	 * */
	Array2Xi boundaryEdgeNodeTags (2, nBoundaryEdges);
	/* now store those edges as node tags,
	 * and store the edge indices associated with each boundary patch. */
	std::vector<ArrayXi> edgeIndicesPerPatch ( nBoundaryPatches );
	
	/* not parallelisable like this, due to dependency on Index "start" */
	Index start {0};
	for ( Index i{0}; i<nBoundaryPatches ; ++i ){
		edgeIndicesPerPatch[i].resize( nBoundaryEdgesPerPatch[i] );
		Index nIndices {0};
		for ( const ElementBlock& eb : msh.elementBlocks ){
			if (
				eb.entityDim == Dim::curve &&
				eb.elementType == ElementType::line &&
				hms::hasEntry( boundaryCurveTags[i], eb.entityTag )
			){
				boundaryEdgeNodeTags.block( 0,start,2, eb.nodeTags.cols() ) =
					eb.nodeTags;
				edgeIndicesPerPatch[i].segment(nIndices, eb.elementTags.size() ) =
					ArrayXi::LinSpaced(
						eb.elementTags.size(),
						start,
						start + eb.elementTags.size() - 1
					);
				nIndices += eb.elementTags.size();
				start += eb.elementTags.size();
			}
		}
	}
	
	return { boundaryEdgeNodeTags, edgeIndicesPerPatch };
}



/**
 * @brief Tags in a MSH file are, if ordered, one-based. They can also be sparse
 * or even arbitrary. The only guarantee is that they are strictly positive.
 * In the hms mesh library, elements are accessed by zero-based consecutive 
 * indices. This function translates tags to indices by looking at the position
 * (i.e. the index) at which a given tag occurs in a list of tags.
 * @param tags the list of tags in which the tag is looked up.
 * @param tag the tag that requires translation.
 * @return the position/index of the passed-in tag within the passed in list
 */
Index tagToIndex( const IdxVector& tags, Index tag ){
	/* if the tags are ordered, they are 1-based, while indices are 0-based, so
	 * a good guess is that the entry is at tag-1. Tags are also strictly
	 * positive, so tag-1 never results in a negative number */
	if (tags[tag-1] == tag) return tag-1;
	
	/* else, the list has to be searched */
	for (uIndex i{0}; i<tags.size(); ++i)
		if (tags[i] == tag) //return static_cast<Index>(i);
			return i;
	
	return hms::signalCode::not_found;
}
/**
 * @overload
 */
Index tagToIndex( Index tag, const IdxVector& tags ){
	return tagToIndex( tags, tag );
}

/**
 * @overload
 */
template<Index rows>
void tagToIndex( Ref<ArrayNX<Index, rows>> m, const IdxVector& tags ){
	OMP_PRAGMA(parallel for schedule(static))
	for ( Index& i : m.reshaped() )
		i = tagToIndex(tags, i);
}

/**
 * @overload
 */
void tagToIndex( IdxVector& v, const IdxVector& tags ){
	OMP_PRAGMA(parallel for schedule(static))
	for ( Index& i : v )
		i = tagToIndex(tags, i);
}

/**
 * @overload
 */
void tagToIndex(
	Ref<Array4Xi> cellNodes,
	const Ref<const ArrayX<char>>& nCellNodes,
	const IdxVector& nodeTags
){
	OMP_PRAGMA(parallel for schedule(static))
	for ( Index col=0; col<cellNodes.cols(); ++col ){
		for ( Index row{0}; row<nCellNodes[col]; ++row ){
			cellNodes(row, col) = tagToIndex( cellNodes(row, col), nodeTags );
		}
	}
}



/**
 * @brief The MSH file format only stores nodes for cells, not edges. However,
 * the fvm method requires edges. This function creates all the inner edges from
 * the existing cell data and records to which cell the new edges belong.
 * @param msh the representation of the MSH file
 * @param domainCellVertices a 4 x n matrix containing the ordered vertex 
 * indices of each domain cell. The MSH file provides ordered node tags, so
 * those merely have to be translated into indices.
 * @param nDomainCellVertices holds the number of vertices (=number of edges)
 * for each domain cell. Used to differentiate between triangles and 
 * quadrilaterals.
 * @param boundaryEdges holds the vertex indices for all boundary edges, i.e. 
 * the translated node tags.
 * @return a pair containing 1. all edges, represented as a matrix of column
 * vectors holding the constituting vertex indices, 2. a matrix holding the edge
 * indices of each cell
 */
auto createDomainCellEdges(
	const MSH& msh,
	const Ref<const Array4Xi>&  domainCellVertices,
	const Ref<const VectorX <char>>&  nDomainCellVertices,
	const Ref<const Array2Xi>& boundaryEdges
) -> std::pair<
	Array2Xi, /* vertex indices of all edges */
	Array4Xi  /* edge indices of domain cells */
>{
/* the number of boundary edges is clear from the MSH file, the number of
 * inner edges can be derived from the number of triangular and quadrilateral
 * cells in conjunction with the number of boundary edges, and the total 
 * number of edges is simply the sum of inner and boundary edges */
	Index nBoundaryEdges { boundaryEdges.cols() };
	Index nInnerEdges { gmsh::nInnerEdges( msh, nBoundaryEdges ) };
	Index nEdges { nInnerEdges + nBoundaryEdges };
	Index nDomainCells { domainCellVertices.cols() };
/* approach for inner edges:
 * somewhat analogous to Mesh::edgesToCells(), a vector of edges at each vertex
 * is built to prevent duplicates. However, here, not all edges are known yet.
 * So this vector is built up in the process of iterating over cells, generating
 * edges and storing parent cells. */
	std::vector<IdxVector> edgesAtNodes ( msh.nNodes );
	
	/* the already existing boundary edge indices are stored at the indices of 
	 * their constituting vertices */
	for ( Index i{0}; i<nBoundaryEdges; ++i ){
		edgesAtNodes[boundaryEdges(0,i)].push_back( nInnerEdges + i );
		edgesAtNodes[boundaryEdges(1,i)].push_back( nInnerEdges + i );
	}
	
	/* consolidate vertex indices of all edges in one variable */
	Array2Xi edgeVertices (2, nEdges);
	edgeVertices.rightCols(nBoundaryEdges) = boundaryEdges;
	
	
	Array4Xi domainCellEdges (4, nDomainCells);
	
	Index edgeIndex {0};
	for ( Index i{0}; i<nDomainCells; ++i ){
		for ( Index j{0}; j<nDomainCellVertices[i]; ++j ){
			/* assume a new edge (block expression doesn't work because of
			 * wrapping around) */
			Vector2i edge {
				domainCellVertices.col(i)[ j],
				domainCellVertices.col(i)[(j+1) % nDomainCellVertices[i] ]
			};
			/* now check whether that edge already exists, otherwise
			 * store it */
			bool edgeIsNew {true};
			for ( Index e : edgesAtNodes[ edge[0] ] ){
				if ( hms::sameEntries(edgeVertices.col(e), edge) ){
					domainCellEdges(j, i) = e;
					edgeIsNew = false;
					break;
				}
			}
			if (edgeIsNew){
				edgeVertices.col(edgeIndex) = edge;
				edgesAtNodes[ edge[0] ].push_back(edgeIndex);
				edgesAtNodes[ edge[1] ].push_back(edgeIndex);
				domainCellEdges(j, i) = edgeIndex;
				++edgeIndex;
			}
		}
	}
	
	return { edgeVertices, domainCellEdges };
}

	


}